Simulations  of  Viscous  Detonations  with 
Detailed  Kinetics  Using  Manifold  and 
Wavelet  Techniques 

by 

Joseph  M.  Powers 
Associate  Professor 

Department  of  Aerospace  and  Mechanical  Engineering 

University  of  Notre  Dame 

presented  at 

Los  Alamos  National  Laboratory 
15  July  1999 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

15  JUL  1999  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1999  to  00-00-1999 

4.  TITLE  AND  SUBTITLE 

Simulations  of  Viscous  Detonations  with  Detailed  Kinetics  Using 

Manifold  and  Wavelet  Techniques 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Notre  Dame, Department  of  Aerospace  and  Mechanical 
Engineering, Notre  Dame, IN, 46556 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

_ _ _  ABSTRACT 

18.  NUMBER  19a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  OS 

unclassified  unclassified  unclassified  Report  (SAR) 

45 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Support 


National  Science  Foundation 
Chemical  and  Transport  Systems,  CTS-9705150, 
Air  Force  Office  of  Scientific  Research, 
Directorate  of  Mathematics  and  Geosciences 
AFOSR  Grant  No.  F49620-98-1-0206, 

and 

Los  Alamos  National  Laboratory 


Acknowledgments 


Prof.  Samuel  Paolucci,  ND-AME 
Dr.  Christopher  Bowman,  Post-Doc,  ND-AME, 
Mr.  Sandeep  Singh,  Ph.D.  Candidate,  ND-AME 
Mr.  Yevgenii  Rastigejev,  Ph.D.  Candidate,  ND-AME 


Outline 


Motivation 

Goals 

Description  of  ILDM  technique 

Summary  of  wavelet  technique 

Detailed  results  for  H2  —  02  detonation 

Strategy  for  HMX  combustion  and  preliminary  results 


Summary 


Motivation 


•  Detailed  finite  rate  kinetics  critical  in  reactive  fluid  mechanics: 

—  Candle  flames, 

—  Atmospheric  chemistry, 

—  Internal  combustion  engines, 

—  Gas  phase  reactions  in  energetic  solid  combustion. 

•  Common  detailed  kinetic  models  are  computationally  expensive. 

—  150  hr  supercomputer  time  for  calculation  of  steady,  laminar, 
axisymmetric,  methane-air  diffusion  flame  (Smooke) 

—  Expense  increases  with 

*  number  of  species  and  reactions  modeled  (linear  effect), 

*  stiffness- ratio  of  slow  to  fast  time  scales,  (geometric  effect). 

—  Fluid  mechanics  time  scales:  10-5  s  to  101  s. 

—  Reaction  time  scales:  ICE14  s  to  102  s. 

•  Reduced  kinetics  necessary  given  current  computational  resources. 

•  Adaptive  discretization  necessary  for  fine  spatial  structures. 

•  Inclusion  of  physical  diffusion  necessary  for  numerical  conver¬ 


gence. 


Why  Diffusion? 

Diffusion  traditionally  not  modelled  in  detonation  studies, 

Argued  that  very  thin  shock  structures,  thickness  =  O(fim),  will 
have  minimal  influence  on  reaction  events, 

However,  inviscid  solutions  to  two-dimensional  reactive  Euler  equa¬ 
tions  in  mildly  unstable  regimes  do  not  appear  to  converge,  while 
viscous  counterparts  do  (Singh,  Powers,  Paolucci,  AIAA-99-0966, 
1999), 

Hypothesis:  inherent  numerical  diffusion  is  selecting  structures  in 
“inviscid”  calculations;  these  evolve  unphysically  with  grid  size, 

When  physical  diffusion  zones  are  resolved  numerically,  grid- 
independent  physical  diffusion  dominates  over  numerical  diffu¬ 
sion. 

Prohibitively  expensive  to  compute  simultaneous  viscous  and  re¬ 
action  zone  structures  with  common  numerical  techniques  and 
actual  physical  parametric  values. 

SPP  modelled  systems  with  reaction  length/diffusion  length  ~  10 
to  achieve  resolved  results;  much  larger  ratios  necessary  to  model 
real  systems. 


Goals 


Implement  robust  new  reduced  kinetic  method  of 

—  Maas,  U.,  and  Pope,  S.  B.,  1992,  “Simplifying  Chemical  Ki¬ 
netics:  Intrinsic  Low- Dimensional  Manifolds  in  Composition 
Space,”  Combust.  Flame ,  88:  239-264. 

—  Lam,  S.  H.,  1993,  “Using  CSP  to  Understand  Complex  Chem¬ 
ical  Kinetics,”  Combust.  Sci.  Tech.,  89:  375-404. 

Extend  method  to  systems  with  time  and  space  dependency. 

Extend  method  to  systems  in  which  fluid  and  chemical  phenom¬ 
ena  evolve  over  similar  time  scales. 

Couple  method  with  new  wavelet  collocation  technique  (Paolucci 
V  Vasilyev)  for  spatial  discretization. 

Applications: 

—  ignition  delay  in  shock  tubes;  detailed  results , 

—  unstable  viscous  detonations, 

—  Bunsen  burner  flames, 

—  rocket  nozzle  flows, 

—  HMX  gas  phase  reactions;  preliminary  manifolds. 


Common  Reduced  Kinetics  Strategies 


•  Fully  frozen  limit:  no  reaction  allowed,  uninteresting 

•  Fully  equilibrated  limit:  commonly  used  in  some  problems 

—  has  value  for  events  in  which  fluid  time  scales  are  slow  with 
respect  to  reaction  time  scales, 

—  misses  events  which  happen  on  chemical  time  scales. 

•  Simple  one  and  two  step  models 

—  require  significant  intuition  and  curve  fitting, 

—  can  give  good  first  order  results, 

—  are  often  not  robust. 

•  Partial  equilibrium  and  steady-state  assumptions 

—  again  require  intuition, 

—  are  not  robust. 

•  Sensitivity  analysis 

—  can  remove  need  to  include  unimportant  reactions, 

—  not  guaranteed  to  remove  stiffness. 


Intrinsic  Low-Dimensional  Manifold  Method  (ILDM) 

•  Uses  a  dynamical  systems  approach, 

•  Does  not  require  imposition  of  ad  hoc  partial  equilibrium  or 
steady  state  assumptions, 

•  Fast  time  scale  phenomena  are  systematically  equilibrated, 

•  Slow  time  scale  phenomena  are  resolved  in  time, 

•  n-species  gives  rise  to  a  n-dimensional  phase  space  (same  as  com¬ 
position  space)  for  isochoric,  isothermal  combustion  in  well  stirred 
reactors, 

•  Identifies  m-dimensional  subspaces  (manifolds),  m  <  n,  embed¬ 
ded  within  the  n-dimensional  phase  space  on  which  slow  time 
scale  events  evolve, 

—  Fast  time  scale  events  rapidly  move  to  the  manifold, 

—  Slow  time  scale  events  move  on  the  manifold. 

•  Computation  time  reduced  by  factor  of  ~  10  for  non-trivial  com¬ 
bustion  problems;  manifold  gives  much  better  roadmap  to  find 
solution  relative  to  general  implicit  solution  techniques  (Norris, 
1998) 


Simplest  Example 


dx 

dt 


=  — lOx, 


dy 

dt 


=  - y , 


Stable  equilibrium  at  (x,y)  = 
ILDM  is  x  =  0 


x(O)  =  x0l 
2/(0)  =  Vo- 

(0,0);  stiffness  ratio  =  10. 


y 


Parameterization  of  manifold:  x(s)  =  0;  y(s)  =  s. 

dy  dy  ds 


■2/(s)  = 


dt  ds  dt ' 
dy  ds 


chain  rule 


ds  dt  ’ 


n\ds 

S  =  (1)di’ 


substitute  from  ODE  and  manifold 
no  longer  stiff! 


5  =  s0e  \ 


x(t)  =  0;  y(t)  =  s0e 


Projection  onto  manifold  for  s0,  induces  small  phase  error. 


Formulation  of  General  Manifolds 


•  A  well  stirred  chemically  reactive  system  is  modeled  by  a  set  of 
non-linear  ordinary  differential  equations: 


dx 

dt 


=  F(x),  x(0)  =  x0, 


x  :  species  concentration;  x  G  3ft 


n 


•  Equilibrium  points  defined  by 


x  =  xeg  such  that  F(xeg)  =  0. 


•  Consider  a  system  near  equilibrium  (the  argument  can  and  must 
be  extended  for  systems  away  from  equilibrium)  with  x  =  x— xeg. 


•  Linearization  gives 


hx 

dt 


=  Fx  x, 


where  Fv  is  a  constant  Jacobian  matrix. 


•  Sclrur  decompose  the  Jacobian  matrix: 


Fx  =  Q  U  Qr 
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Formulation  of  General  Manifolds  (cont.) 


Q  is  an  orthogonal  matrix  with  real  Schnr  vectors  q,  in  its  columns. 

U  is  an  upper  triangular  matrix  with  eigenvalues  of  Fx  on  its 
diagonal,  sometimes  placed  in  order  of  decreasing  magnitude. 

The  Schnr  vectors  qt  form  an  orthonormal  basis  which  spans  the 
phase  space, 

We  then  define  m  slow  time  scales,  m  <  n. 


Next  define  a  non-square  matrix  W  which  has  in  its  rows  the 


Schur  vectors  associated  with  the  fast  time  scales: 

( .  ub  ,  i  . \ 


W 


7 

Qm+ 1 
Qm+ 2 


\  •  •  •  •  •  •  n ^  ...  ...  I 

\  Hn  / 

Letting  the  fast  time  scale  events  equilibrate  defines  the  manifold: 


W  •  F(x)  =  0. 


If  m  =  0,  no  slow  time  scales,  W  =  Qr,  and  W  •  F(x)  =  0 
implies  Qr  •  F(x)  =  0,  implies  F(x)  =  0:  the  equilibrium  point 
is  the  low  dimensional  manifold! 


A  Simple  Example 


Consider 


dx 

dt 

dy_ 

dt 

Equilibrium  points: 

/O' 


—  lOOx  +  y  sin?/,  x(0)  =  x 

x3  -  y,  2/(0)  =  2/0- 


F  = 


(  — lOOx  +  y  siny 


V 


x3  —  y 


1  x 

\y. 


(o' 
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Other  equilibrium  points  exist! 


Near  the  equilibrium  point  (0,0),  linearization  gives 

(  x 


/  dx\ 

t 

dt 

rfy  , 
\  dt  / 

V 

0 


-1 


\y. 


which  is  obviously  stable. 
Schur  decomposition  is  trivial: 


/-100  0 

V  0  -1 

Form  the  manifold: 
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The  ILDM! 


A  Simple  Example 


—  =  -100x  +  ysin(y)  x(0)  =  xo 
^=x3-y  y(0)  =  yo 


Simple  Example:  Parameterization  and  Stiffness  Reduction 


=  —  lOO.r  +  ?/ sin  2/ , 


—  =  x  —y. 
dt 


•  Time  scales  near  origin:  t\  =  1.0,  =  0.01.  Stiff. 

•  First  approximation  to  manifold  is  x  =  -j^ysiny. 


•  Parameterize  manifold  as 


•  Chain  rule  gives 


x  =  ioPsins’ 


y  =  s. 

dy  dy  ds 
dt  ds  dt 


•  Substitute  from  ODEs  and  parameterization: 

3  dy(s)  ds 

X  (S) 

— -s6  sim  s  —  s  =  (1) — , 
106  (It7 

ds  loo 
—  =  sin  s  —  s 
dt  106 


•  Linearize  near  equilibrium  at  origin: 


Time  scale:  r  =  1.0  No  longer  stiff! 


Solve  ODE  for  s(t),  substitute  to  get  x(s(t)),  y(s(t)): 


x  ~ - sQe  sin 

-1  u 


in  (s0e  r) 


,  y  ~  s0e 


Example:  Zeldovich  Mechanism  of  NO  Formation 

•  Mechanism  (two  elements,  five  species,  two  reactions): 


1.  0  +  N2  -►  NO  +  N,  =  1.8  x  1014  ^  exp  p383r70  A'), 

2.  NO  +  N  ->  O  +  N2,  k2  =  3.8  x  1013  ^exp(^f-^), 

3.  N+02  -  NO+O,  k3  =  l.SxlO10  ^-Texp 

4.  NO+O  -  ]V+02,  k4  =  3.8xW9^1?Texp  (~208r20  *), 

•  Take  T  =  1400  K,  then 

1.  fci  =  2.252  x  102 

1  mol  S 

2.  fc2  =  2.805  x  1013 

3.  fc3  =  8.905  x  1011 

d  mol  s 

4.  k4  =  1.851  x  106 

^  mol  S 

•  Law  of  mass  action  for  [A^] ,  for  example,  gives 

®  =  —  *i[JV2][0]  +  k2[NO][N}. 

(If 

•  For  all  species,  law  of  mass  action  yields  five  non-linear  ODEs: 

h[N2}[0} 
k2[N][NO] 
h[N}[02] 
h[NO][0] 


Example:  Zeldovich  Mechanism  of  NO  Formation,  cont. 


•  To  elucidate  naturally  conserved  variables,  use  elementary  row 
operations  to  cast  system  in  non-unique  row  echelon  form: 


d 

dt 


'  in 

\ 

(1 

-1 

-1 

1  ^ 

[NO]  -  [iV] 

0 

0 

2 

-2 
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0 
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0 
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0 

0 
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[N], 

v0 

0 

0 

0  y 

(  h[N2}[0]  \ 
k2[N][NO 
h[N}{02] 
U4  [N0\[0\j 


•  We  are  left  with 


—  two  ODEs 

—  three  algebraic  constraints:  conservation  of  N  atoms,  O  atoms, 
and  number  of  molecules 

—  easily  reduced  to  two  ODEs  in  two  unknowns:  [TV],  [NO]. 

•  We  will  reduce  the  two  ODEs  to  one  ODE  by  imposing  the  man¬ 
ifold  equation  W  •  F(x)  =  0,  effectively  equilibrating  the  fast 
time  scale. 


Example:  Zeldovich  Mechanism  of  NO  Formation,  cont. 


•  Consider  first  the  intrinsic  algebraic  constraints: 

/ 2[tV2]  +  [NO]  +  [N]  \  /0\ 

I  +  =  0  • 

Uo2]  +  [iVO]-[JV]j  Vo/ 

•  Integrate  these  equations: 

2[JV2]  +  [WO]  +  [N]  =  Ci, 

[o]  +  [JV]  =  c2, 

2[02]  +  [NO]  -  [N]  =  C3. 

The  constants  Ci,  C2,  C3  come  from  initial  conditions. 

•  Solve  equations  for  secondary  variables  in  terms  of  [N] ,  [NO] : 

[AVI  =  \  {Cl  -  [ATO]  -  [N]) 

[O]  =  C2-[N] 

[02]  =  \  (C3  -  [NO]  +  [AT]) 

•  Note  that  rearrangement  of  the  algebraic  constraints  demonstrates  ele¬ 
ment  and  molecule  conservation: 

2[N2]  +  [N]  +  [NO]  =  Ch 
2[02]  +  [N0\  +  [0]  =  C2  +  C3, 

[iV]  +  [ATO]  +  [Ay  +  [O]  +  [02]  =  Cl  +  C'3+  c2. 


Example:  Zeldovich  Mechanism  of  NO  Formation,  cont. 


Substitution  of  algebraic  constraints  into  ODEs  for  [. N ]  and  [NO]  gives 
two  autonomous  ODEs  well-suited  for  dynamic  systems  analysis: 

cM  =  ^ct-m^-m-iNO]) 

-k2[N][NO] 

~W\  (C3  +  [iV]  -  [JVO]) 

+ki[NO]  (C2  -  [ N ]) 

CMA  =  |  (c?2  -  [JV])  (Cl  -  [AT]  -  [JVO]) 

-fc2[Af][AfO] 

+yM  (C3  +  [iV]  -  [1VO]) 

-kt[NO}  (C2  -  [AT]) 


Take  as  initial  conditions 

777  n/p 

[fV]  =  [fVO]  =  [N2\  =  [O]  =  [02]  =  0.001 - -. 


Equilibrium  when  right  hand  side  zero 


Three  roots-one  physical,  two  unphysical: 


m 

\[NO] 


.16  x  10~n  ^§-\  /-1.15  x  10"11  (  —2.00  x  10 

CTTI  *  r-rn J  1  I 


>.78  x  10"6  ^ 


-2.78  x  10~6  ^ 

cm 6 


^—3  mole 


cm- 


0.00  x  10°  ^ 

cm0 


Example:  Zeldovich  Mechanism  of  NO  Formation,  cont. 


•  Linearization  of  equations  near  physical  equilibrium  gives 
d  (  [N}~  1.16  x  1CT11  \  /  —9.67  x  108  3.38  x  1( 


[NO]  -  2.78  x  10 


d  (  [N]  -  1.16  x  1CT11 
dt  V  [iVO]  -2.78  x  10”6 


-9.67  x  10s  3.38  x  106 

8.11  x  108  -4.03  x  103 

[IV]  -  1.16  x  10~n  \ 

[NO]  -  2.78  x  10"6/ 
-0.766  -0.643  \ 

Q 

0.643  -0.766 ) 

-9.67  x  108  3.38  x  103 

0  -1.19  x  103 

-0.766  0.643  \ 

qt 

-0.643  -0.766 ) 

[IV]  -  1.16  x  10~n  \ 

[NO]  -  2.78  x  lO"6/ 


•  Condition  number  (stiffness  ratio)  =  =  8.1  x  105. 


•  Locally  the  ILDM  is  defined  by 


-0.766  0.643 


W-F(x)  =  0, 

n([Ar].[ivo])\ 


F2([iV],[iVO]) 


=  0. 


•  Use  arc  length  continuation  methods  to  define  complete  ILDM 

•  The  physical  equilibrium  has  negative  eigenvalues:  stable. 


•  The  non-physical  equilibria  have  positive  eigenvalues:  unstable. 


10' 


Adaptive  Multilevel  Wavelet  Collocation  Technique 


•  Summary  of  standard  spatial  discretization  techniques 

—  Finite  difference-good  spatial  localization,  poor  spectral  local¬ 
ization,  and  slow  convergence, 

—  Finite  element-  good  spatial  localization,  poor  spectral  local¬ 
ization,  and  slow  convergence, 

—  Spectral-good  spectral  localization,  poor  spatial  localization, 
but  fast  convergence. 

•  Wavelet  technique 

—  See  e.g.  Vasilyev  and  Paolncci,  “A  Fast  Adaptive  Wavelet  Col¬ 
location  Algorithm  for  Multidimensional  PDEs,”  J.  Comp. 
Phys .,  1997, 

—  Basis  functions  have  compact  support, 

—  Good  spatial  localization,  good  spectral  localization,  and  fast 
convergence, 

—  Easily  formulated  to  adapt  spatially  to  capture  steep  gradients 
via  adding  collocation  points, 

—  Spatial  adaptation  is  automatically  and  dynamically  adaptive 
to  achieve  prescribed  error  tolerance. 


Ignition  Delay  in  Premixed  H2-O2 


Consider  standard  problem  of  Fedkiw,  Merriman,  and  Osher,  J. 
Comp.  Phys .,  1996, 

Shock  tube  with  premixed  iO,  0%  and  Ar  in  2/1/7  molar  ratio, 
Initial  inert  shock  propagating  in  tube, 

Reaction  commences  shortly  after  reflection  off  end  wall, 
Detonation  soon  develops, 

Model  assumptions 

—  One-dimensional, 

—  No  diffusion  (one  case);  mass,  momentum,  and  energy  diffu¬ 
sion  (another  case), 

—  Nine  species,  thirty-seven  reactions, 

—  Ideal  gases  with  variable  specific  heats. 


Compressible  Reactive  Navier-Stokes  Equations  for  H2-O2  Problem 


dp  d  ,  N 

at  +  m (pu)  =  °’ 


mass 


U  .  .  1/  /  o  n  \ 

—  (pw)  +  —  [pu  +  P  —  t  =  0,  momentum 
at  ax  v  ' 


5 

dt 


d_ 

dx 


+  u(P 


energy 


d  d  M  f—E-\  N  ( oYi\Vk3 

n(j^)  ■  sPecies 


n  y 

P  =  p^RT  V  — thermal  equation  of  state 

ti 


^  /  rT  ~  P 

e  =  ^2  Yt  +  Jt  cpi(T)dT J - ,  caloric  equation  of  state 


4  5m 

r  = 


Newtonian  gas  with  Stokes’  assumption 


N  QY 

=  -p  £  A  J 

3= 1 


5x 


Fick’s  law 


9  = 


augmented  Fourier’s  law. 


N  =  9  species:  FT2,  02,  H,  O,  OH,  H202,  H20,  H02 ,  Nr 
M  =  37  reactions 


Operator  Splitting  Technique 


•  Equations  are  of  form 


d  d 

fa  qOM)  +  ^f(qOO))  =  g(qOM))- 

where 

q  =  (d,  pu,  p  +  y  j  ,  pYi 

•  f  models  convection  and  diffusion 

•  g  models  reaction  source  terms 

•  Splitting 

1.  Inert  convection  diffusion  step: 

+  J^f(qOM))  =  o, 

^q  i(t)  =  -Aa;f(q  i{t)). 

Ax  is  either  Godunov  or  wavelet  discretization  operator. 

2.  Reaction  source  term  step: 


d 

— q(x,t)  =  g(q(>G)), 
^q*W  =  g(q*W)- 


•  Operator  splitting  with  implicit  stiff  source  solution  can  induce  non¬ 
physical  wave  speeds!  (LeVeque  and  Yee,  JCP  1990) 


ILDM  Implementation  in  Operator  Splitting 


•  Form  of  equations  in  source  term  step: 


p  \ 

(0\ 

d 

pu 

0 

dt 

P(e  +  i) 

0 

\  pYi  ) 

•  Equations  reduce  to 

P  Poi  ^  Uo,  6 

dYi  uj 
dt  pa 

•  uj  has  dependency  on  p,  e,  and  Y% 

•  ODEs  for  Y{  can  be  attacked  with  manifold  methods  when  man¬ 
ifold  with  p,  e,  H  and  O  parameterization  is  available. 

•  In  premixed  problem,  H  and  O  element  concentrations  are  re¬ 
markably  constant,  reducing  the  dimension  by  two! 

•  Full  equations  integrated  until  sufficiently  close  to  manifold 

•  Once  on  manifold,  simple  projection  used  to  return  to  manifold 
following  convection-diffusion  step 


Sample  ILDM  for  H2  -  02 


Projection  of  ILDM  in  P20,  H202  plane, 

Adiabatic  (e  =  525  kJ/kg ),  isochoric  (p  =  0.25  %/m3),  element 
concentrations  of  H  and  O  constant, 

Complete  manifold  tabulated  in  three  dimensions:  p,  e,  Yh2o-> 

So  we  have  e.g.  P  (p,  e,  lp2o) ,  T  (p,  e,  lp2o) ,  YH  (p,  e,  lp2o) , . . . 
Linear  interpolation  used  for  points  not  in  table, 

Captures  ^0.1  ps  reaction  events. 


H20  mass  fraction 


Inviscid  H2  —  O2  Ignition  Delay  with  and  without  ILDM 


•  No  diffusion, 

•  Godunov  spatial  discretization,  400  uniform  finite  difference  cells, 

•  Implicit  (trapezoidal)  convection  step;  Implicit  (dlsode)  or  ILDM 
reaction  step, 

•  Correction  of  Fedkiw  adopted  to  suppress  artificial  entropy  layer 
after  shock  reflection  (see  Menikoff,  1994). 
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Viscous  H)  —  O2  Ignition  Delay  with  Wavelets 

•  Mass,  momentum,  and  energy  diffusion  modelled, 

•  Wavelet  spatial  discretization,  explicit  convection-diffusion  time  step¬ 
ping,  implicit  reaction  time  stepping, 

•  300  collocation  points,  15  wavelet  levels, 

•  Viscous  shocks,  induction  zones,  and  entropy  layers  spatially  resolved! 

•  t  =  180  /is. 


remperature  (K) 


Viscous  H2  —  O2 


•  t  =  190  /is 


x  105 


Density  (kg/m  ) 


1  Delay  with  Wavelets 


remperature  (K) 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  =  200  fis 


x  105 


remperature  (K) 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  =  230  /is 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  —  180  fis 

•  species  mass  fractions  plotted  vs.  distance 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  =  190  fis 

•  species  mass  fractions  plotted  vs.  distance 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  =  200  fis 

•  species  mass  fractions  plotted  vs.  distance 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets 


•  t  =  230  /is 

•  species  mass  fractions  plotted  vs.  distance 


Temperature  (K) 


Post  Reflection  Entropy  Layer?:  Viscous  Wavelet  Results 

•  No  significant  entropy  layer  evident  on  macroscale  after  shock 
reflection  when  resolved  viscous  terms  considered, 

•  Inviscid  codes  with  coarse  gridding  introduce  a  larger  entropy 
layer  due  to  numerical  diffusion, 

•  Unless  suppressed,  unphysically  accelerates  reaction  rate. 


Pressure  (Pa)  Temperature  (K) 


Post  Reflection  Entropy  Layer:  Viscous  Wavelet  Results 


•  small  entropy  layer  evident  on  finer  scale, 

•  temperature  rise  ~  5  K\  dissipates  quickly, 

•  inviscid  calculations  before  adjustment  give  persistent  tempera¬ 
ture  rise  of  -  20  K\  reaction  acceleration  small. 
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Viscous  H_2  —  O2  Ignition  Delay  with  Wavelets 
Close-up:  Viscous  Shock  Stucture  and  Induction  Zone 

•  t  =  230  (is , 

•  Induction  zone  length:  ^  470  nmi 

•  No  significant  reaction  in  viscous  shock  zone. 
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Viscous  H_2  —  O2  Ignition  Delay  with  Wavelets 
Closer-up:  Viscous  Shock  Stucture  Only 

•  t  =  230  (is 

•  predicted  shock  thickness:  ~  50  /im. 


x  105 


Level  Level 


Viscous  H2  —  O2  Ignition  Delay  with  Wavelets, 
Instantaneous  Distributions  of  Collocation  Points 

•  t  =  180  fis,  two-shock  structure  with  consequent  collocation 
point  distribution, 

•  t  =  230  [is,  one-shock  structure  with  evolved  collocation  point 
distribution. 
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Application  to  Gas  Phase  HMX  System 

Simulating  isobaric  HMX  combustion  computationally  intensive, 

Most  effort  in  solving  gas  phase  convection,  reaction,  diffusion, 

Based  on  45  species,  232  step  mechanism  of  Yetter,  et  ah, 

Fastest  time  scales  predicted  ICE16  s  (non-physical?), 

Stiffness  ratio  1011  (vs.  109  for  H2  —  O2), 

Equations  for  gas  phase  combustion  of  HMX  are  of  form 

d  d 

^q(M)  +  g^f(q  OM))  =  g(q(M)), 

Adiabatic,  isobaric, 

Operator  splitting  appropriate, 

For  non-premixed  problem,  higher  dimension  (>  8?!)  manifolds 
necessary! 

Will  need  to  parameterize  by  (/i,  p,  H,  O,  TV,  C,  Mr,  >  one  free  parameter) 
107  <  h  <  1011  erg/ g\  1CT5  <  p  <  1CT3  p/cm3;  10-32  <  XAr  <  1CT2; 

0  <  Xc  <  101;  0  <  Xh  <  101;  0  <  xn  <  101;  0  <  Xo  <  101. 

(Liau,  1999) 

Three-dimensional  manifold  for  preliminary  premixed  problem? 


ILDM  for  Gas  Phase  HMX  System 


Based  on  45  species,  232  step  mechanism  of  Yetter,  et  ah, 
Adiabatic  (h  =  62  x  109  erg/ g)  ,  isobaric  (P  =  32  bar), 
projection  in  Yjv2,  Yc02  plane- 


Summary 


•  Robust  method  in  place  to  compute  manifolds  with  arbitrary 
variables  held  constant  (e.g.  P,  p,  h), 

•  Effort  still  needed  on  improving  technique  of  projecting  onto  man¬ 
ifold  initially, 

•  Fast  linear  interpolation  scheme  in  place  for  table  lookup, 

•  Robust  method  in  place  to  solve  less  stiff  differential  equations 
on  or  near  manifold, 

•  Operator  splitting  allows  implementation  of  manifold  in  solving 
PDEs, 

•  Adaptive  multilevel  wavelet  collocation  method  gives  dramatic 
spatial  resolution, 

•  Full  coupling  of  ILDM  and  wavelet  methods  soon  forthcoming, 

•  Detailed  studies  of  efficiency  improvement  necessary, 

•  More  general  manifold  techniques  need  developed  to  allow  strong 
fluid-chemistry  coupling  and  relaxation  of  eigenmodes  to  steady 
state  solutions. 


